Effect of the cluster size in modeling the H2 desorption and dissociative adsorption on 

Si(OOl) 
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Three different clusters, SigHi2, SiisHie, and Si2iH2o, are used in density-functional theory calcu- 
lations in conjunction with ab initio pseudopotentials to study how the energetics of H2 dissociative 
adsorption on and associative desorption from Si(OOl) depends on the cluster size. The results are 
compared to five-layer slab calculations using the same pseudopotentials and high quality plane- 
wave basis set. Several exchange-correlation functionals are employed. Our analysis suggests that 
the smaller clusters generally overestimate the activation barriers and reaction energy. The Si2iH2o 
cluster, however, is found to predict reaction energetics, with Ea CB = 56±3 kcal/mol (2.4±0.1 eV), 
reasonably close (though still different) to that obtained from the slab calculations. Differences in 
the calculated activation energies are discussed in relation to the efficiency of clusters to describe 
the properties of the clean Si(001)-2xl surface. 

PACS numbers: 68.35.-p, 82.65.My 
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I. INTRODUCTION 

Cluster models are a frequently used tool for study- 
ing different aspects of physics and chemistry of clean 
surfaces, adsorption, and surface chemical reactions. An 
example which attracted considerable experimental and 
theoretical interest over the last years is the H dissocia- 
tive adsorptkm and associative desorption on the Si(OOl)- 
2x1 surfaced which is a subject of this study. The in- 
triguing experimental result that .t|h|G H2 desorption from 
Si(OOl) follows first-order kineticsata has triggered an in- 
tense theoretical activity in this field mainly concentrated 
on the mechanism(s) leading to such an unusual behav- 
ior. The available first principles calculations address 
the latter question on the basis of two different models: 
(i) the cluster approximation!, using either configuration- 
interaction_(CI) methodsETt3 or density-functional the- 
ory (DFT)Bt3 to describe the exchange and correlation 
effects or (ii) extender! slab models for the Si(001)-2xl 
surface using DFThXH. 

The DFT slab calculations all agree in their conclu- 
sions supporting the pre-pairing mechanismcl according 
to which two hydrogens are pre-paired on the same Si 
surface dimer and associatively desorb through an asym- 
metric transition state (TS). The cluster calculations, 
however, have led to different conclusions. All these cal- 
culations find a rather high barrier for desorption of two 
hydrogens from a -single Si dimer, e.g. £^ es ._j= 74-75 
kcal/mol (3.2 eV)B, 85,-86 kcal/mol (3.7 eV)E3, 82-85 
kcal/mol (3.6-3.7 eV)E3. Comparing to the experimen- 
tal ftftpftti 011 energy of desorption of ~ 58 kcal/mol (2.5 
eV)0aO'E3, these findings were interpreted as being com- 
pelling evidence against the pre-paring mechanism. In 
an attempt to reconcile the experimentally observed en- 
ergetics and kinetics of desorption from the monxihydride 
phase, various defect-mediated mechanismstTHiij were 
suggested including formation of metastable dihydride 



species as an intermediate step. All the above studies 
have in common that their argumentation rests on com- 
putational schemes based on CI and the small and simple 
SigHi2 cluster to model the ij>i(00T)-2x 1 surface except 
the work by Nachtigall et al.U where DFT has been ^m- 
ployed.-_.In some studies larger clusters, like SiigH^B or 
Si2iHji_, have been used. Nevertheless, the effect of clus- 
ter size on the energetics of H2 adsorption/desorption on 
Si(OOl) has not been analyzed in detail in the literature. 

The only cluster-based support to direct desorption 
via the pre-pairing mechanism-.comes from the DFT 
calculations by Pai and Doren__. Using—the non-local 
Becke-Lee- Yang-Parr (BLYP) functionat_]'__, they find 
E^ es = 64.9 kcal/mol (2.8 eV) including a zero-point en- 
ergy (ZPE) correction. Given the uncertainties originat- 
ing both from the use of different functionals and from 
inherent limitations of the cluster approximation, they 
considered their result to be compatible with the exper- 
imental desorption energy. The first possible source of 
error, the reliability of the functional, was addressed by 
Nachtigall et alLH, who systematically compared various 
density functionals to the results of the state-of-the-art 
computational tool in quantum chemistry, an extrapo- 
lated quadratic CI method. For the sake of computa- 
tional feasibility, they concentrated on a few simple test 
cases, four reactions involving silanes, and a Si2Hg clus- 
ter with a geometry chosen to mimic H2 desorption from 
Si(OOl). While the extrapolated quadratic CI method 
gives a reference value of 90.4 kcal/mol (£4)2 eV) for 
E^ cs , the Pcrdew-Wang (PW91) functional^! underesti- 
mates E^ cs by 9.5 kcal/mol (0.41 eV) compaped-to this 
reference. The Becke-Perdew (BP) functional-]'-- gives 
an even lower barrier, E^ os = 79.4 kcal/mol (3.44 e\CL 
Generally, Nachtigall et al. find the B3LYP functional-- 
to give closest agreement with their CI calculations, while 
the BLYP functional is performing second best. A similar 
trend concerning the performance of different function- 
als was also observed in the case of H diffusion on the 
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Si(001)-2xl surfaced. However, we find it not to be a 
priori clear whether their conclusions could be applied 
to better cluster approximations to the Si(OOl) surface, 
because the electronic wave functions at a surface are 
generally more extended than in a cluster and this may 
naturally affect the electron-electron correlations. It is 
one of the aims of the present paper to study the perfor- 
mance of different functionals with the size of the clus- 
ters. 

A second issue which hampers a clear-cut comparison 
of different approaches lies in the geometric structure of 
the clean Si(OOl) surface. While the DFT calculations 
give the correct description of the geometry of the clean 
surface, normally .a a (2 x 2) or c(4 x 2) reconstruction with 
buckled Si dimcrsESEH, CI calculations predict a symmet- 
ric ground state of the clusters. This could imply sub- 
stantial differences between the two approaches in the 
surface relaxation during the adsorption/desorption pro- 
cess as well as in the surface electronic structure. Since 
H 2 dissociatjpn-jm Si(OOl) is known to couple to the sur- 
face motionolijlij, one could expect these differences in 
the description of the clean Si(OOl) surface also to influ- 
ence the reaction energetics predicted by different calcu- 
lational approaches. 

In the light of the above we find of particular interest to 
bridge the DFT slab and cluster calculations by studying 
the convergence of the H2 adsorption/desorption energet- 
ics with cluster size. Since the results of the DFT calcula- 
tions agree about a direct desorption-orocess and the very 
recent experiment by Flowers et alE3 also supports this 
mechanism, we shall perform our analysis adopting the 
pre-pairing scenario for the Si(OOl) surface. A descrip- 
tion of the computational method used in the present 
work is given in the next section. In Sec. Ill particu- 



lar attention is paid to the structure of the clean Si(OOI) 
surface. Clusters with one, two and three surface dimers 
are used in Sec. to study the energetics of the adsorp- 
tion/desorption process. A summary of the results and 
discussion is presented in the last section. 



Each of them represents the topmost four layers of the 
reconstructed surface. They differ by the number of Si- 
Si surface dimers they contain. The larger SiisHig and 
Si2iH2o clusters are derived from SigHi2 by adding one 
and two dimers, respectively, in the [110] direction along 
with their full coordination of second layer Si atoms, one 
third and one fourth layer Si atom. All dangling bonds of 
the silicons in the subsurface layers are saturated with hy- 
drogen atoms. By increasing the cluster size in this way 
we aim at reaching three ascendingly improved approxi- 
mations to the clean Si(001) surface: SigHi2 is the mini- 
mal one to represent the symmetric 2x1 reconstruction; 
Sii5Hi 6 is the smallest cluster that enables to model the 
p(2 x 2) surface reconstruction. Finally, Si2iH2o contains 
one surface dimer surrounded by two others, thus having 
the same local environment as on the surface. In the con- 
text of H2 desorption from the monohydride phase, the 
larger clusters also allow to study the interaction between 
adjacent occupied dimers. 



Si,,H 




FIG. 1. BP optimized geometries of (a) the clusters used 
to model the Si(001)-2x 1 surface and (b) the five-layer p(2x 2) 
slab used as a reference system. 



II. CALCULATIONS 



A. The systems used for modeling H2/Si(001) 

Generally, cluster models for surface chemical pro- 
cesses are treated in conjunction with a basis set for the 
electronic states consisting of localized orbitals. How- 
ever, special care is required in choosing a basis set which 
meets the desired level of accuracy. In order to isolate 
the different approximations inherent in the choice of the 
cluster and the basis set size, we decide to perform to- 
tal energy calculations within the.. DFT scheme as im- 
plemented in the fhimd packagd23, employing a plane 
waves basis set with well-controlled convergence prop- 
erties. We report results for the three different clus- 
ters shown in Fig. |l| (a), SigHi2, Sii5Hi6 and Si2iH2o- 



Since the super cell approach implemented in fhimd 
implies periodic boundary conditions for the wave func- 
tions, the clusters are placed in a sufficiently large or- 
thorhombic unit cell, where they are separated by w 7 
A to avoid unwanted interactions. Since the electronic 
states of the cluster Hamiltonian have no dispersion, it 
is sufficient to calculate them at a single point in the 
Brillouin zone (BZ), k r = (0,0,0). 

To link the current discussion with the available slab 
calculations, we perform as a final step a set of calcu- 
lations employing a p(2 x 2) slab with five Si layers, 
Fig. g (b), similar to Ref. The fc-space integration is 
performed using 16 k|| points in the whole surface Bril- 
louin zone of the 2x2 unit cell. 
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Both cluster and slab calculations are carried out using 
the local-density approximation (LDA) to the exchange- 
correlation functional, as obtained form the Monte Carlo 
results by Ceperley and, AlderEil in the parameterization 
of Perdew and Zungerca. For comparison all calculations 
are repeated with gradient corrections as described by the 
BP, PW91 and BLYP functionals. Our choice to include 
these functionals in the test is partially motivated by the 
good agreement between previous slab calculations and 
the experiment, but on the other hand by the shortcom- 
ings of the BP and PW91 functionals, cf. Ref. [ll], found 
in the case of small systems. DFT plane-wave calcula- 
tions employing either slabs or Sig or larger clusters and 
the BLYP functional have not yet been reported for the 
H2/Si(001) reaction energetics. 

We employ norm-conserving sp-nonlocal pseudopoten- 
tials for Si atoms, generated from a separate all-electron 
calculatiop_pf the Si atom for each exchange-cacrelation 
functional^ according to Hamann's schemea. For 
gradient-corrected DFT calculations, the use of these 
consistently constructed pseudopotentials ensures the 
proper description of core-valence exchange within each 
of the gradient-corrected functionals usedlH For the hy- 
drogens passivating Si dangling bonds, a s-nonlocal pseu- 
dopotential is generated following the Troullier and Mar- 
tins prescriptiont3. However, for the hydrogen atoms 
taking part in the reaction, we employ the full 1/r po- 
tential. For all energies quoted in the paper, plane waves 
with a kinetic energy up to E cut — 30 Ry (408 eV) were 
included in the basis set. While geometries and relative 
energies for systems consisting entirely of Si atoms are 
well converged already at 18 Ry (245 eV), the high qual- 
ity basis set is required for a correct description of the 
hydrogen wave functions close to the core. 



B. Structure optimization 



structionEZI. To distinguish elastic from electronic effects, 
we employ a two-step procedure to determine the equi- 
librium structures of the bare clusters. First we sample 
the total energy as a function of the dimer buckling an- 
gle a, as done, for example, in Refs. 
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but without 

relaxation of any of the deeper layers. Thus any elastic 
interactions between Si dimers are avoided. Neighbor- 
ing dimers (if there are any) are buckled in opposite di- 
rections, with buckling angles a and —a, to mimic the 
p(2 x 2) surface reconstruction. 

The structure of the bare clusters is finally determined 
by unconstrained relaxation of the topmost two Si lay- 
ers. Here we use the energy minimum as a function of 
a determined in the previous calculations as input for 
the starting geometry. The two pairs of hydrogens sat- 
urating the second layer Si bonds in the [110] direction, 
which would correspond to adjacent Si dimers on the 
Si(001) surface, are also allowed to relax. We have also 
tested full relaxation of all layers in the case of the SigH 12 
and Sii5Hi 6 clusters within LDA. The lack of any geo- 
metric constraints, however, tends to overestimate the 
surface relaxation and introduces unrealistic atomic dis- 
placements. Full relaxation would only be appropriate 
for clusters studied as objects in their own interest, rather 
than as an approximation to the Si(001) surface. 

The structure of the monohydride phase is determined 
by relaxing the adsorbate and the two Si layers beneath 
it. In the case of the S121 cluster the H2 molecule is 
adsorbed on the middle dimer. The geometries of tran- 
sition states (TS) are determined by a search algorithm 
using the ridge method proposed by Ionova and CarterE3. 
Since geometries are less sensitive to the qualtity of the 
basis set than the total energies, we have found it suffi- 
cient to perform geometry optimizations at a plane- wave 
cut-off of E cut = 18 Ry. The structures are considered 
converged when all forces are smaller than 0.05 eV/A. 



We perform separate structure optimizations for each 
cluster size and for the LDA, BP, PW91 and BLYP func- 
tionals. The first step in all structure relaxation runs was 
to optimize the positions of the terminating hydrogens, 
with the subsurface silicon coordinates kept fixed at their 
bulk values, and the symmetric dimer bond length set to 
d = 2.28 A as found in Ref. ||. 

One of our present goals is to investigate to what extent 
clusters are appropriate for a description of the Si(001) 
surface. One important feature of the clean surface is the 
buckling of the Si surface dimers (see the next Section). 
It is therefore interesting to test if the buckling can be 
modeled by clusters of appropriate size. The failure of 
small clusters to reproduce the buckled surface structure 
is sometimes attributed to the neglect or underestimate 
of elastic interactions between Si dimers in this approach. 
However, the importance of elastic interactions is estab- 
lished only for the alternation of buckling angle ('anti- 
ferromagnetic ordering') in the p(2 x 2) or c(4 x 2) recon- 



III. THE CLEAN Si(001)-2xl SURFACE 

On the Si(001) surface, the surface Si atoms form 
dimers, leading to the (symmetric) 2x1 reconstruction. 
The surface can reduce its symmetry if the Si surface 
atoms relax to different heights, i.e. the dimers are buck- 
led. This leads to the formation of lower symmetry pat- 
terns, the asymmetric p{2 x 1), p{2 x 2), c(4 x 2). These 
reconstructions are characterized by the dimer buckling 
angle a, the dimer bond length d, and the energy favor 
per dimer AE with respect to the symmetric 2x1 re- 
construction. 

The calculations for H2 desorption from Si(001) using 
slabs together with DFT and those based on many-body 
wave functions (like CI) already differ in their description 
of the clean surface geometry. In the latter calculations, 
mostly the SigHi2 cluster with one symmetric Si dimer 
was used as reference state for the clean Si(001) surface. 
Radeke and Carter verified that the symmetric cluster is 
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the geometric-ground state at the level of theacy used 
in their workE3. Recent theoretical studiesOcil on Si 
clusters come to conflicting conclusions and thus leave 
the discussion about the ground state symmetry of Si 
clusters still open. 

Oi^|-the,|— other hand, the DFT slab calcula- 
tionalall3'E§E3e3 are consistent in their predictions and 
show good agreement with recent experiments. Direct 
evidence forj-the buckling comes from low-temperature 
STM imagestZl and structure determinations by LEEEO. 
Furthermore, observed core-level shiftscj are inconsistent 
with symmetric-surface dimers, but can be explained by 
buckled dimersEJ. The buckled surface reconstruction 
has gained further support by the good agreement bp. 
tween the measured dispersion of surface band-statescll 
with calculations using the GW approximationEl 

The different symmetries of reconstruction correspond 
to different electron,ic || st | :pjritures at the surface. Previ- 
ous theoretical workli3oa has established the following 
picture: After dimerization of the surface Si atoms, they 
would both remain to have dangling bonds occupied by 
one electron, i.e. a degenerate electronic ground state. 
There are two principal possibilities for a lowering of the 
electronic energy. Firstly, a bonding and antibonding lin- 
ear combination of orbitals could be formed (similar to a 
7r-bond in a free dimer), only one of which is occupied. 
In this case the Si surface dimer would remain symmet- 
ric. The second possibility is a Jahn- Teller-like splitting 
of the degeneracy. By buckling the Si dimer, the lower Si 
atom comes to an almost planar bonding configuration 
with its three neighbors, while the upper Si atom reaches 
a pyramidal configuration. A rehybridization of the or- 
bitals at each of the Si atoms results in a lowering of the 
dangling orbital at the upper Si atom, and in an up-shift 
of the orbital at the lower Si atom, which is accompanied 
by a transfer of electron density from the lower to the up- 
per atom. The preferred way of stabilization depends on 
several factors: the possible strength of the 7r-bond, the 
possible energy gain due to rehybridization, the ability of 
the system to screen the increased Coulomb repulsion in 
the dangling bond of the upper Si atom, and the energetic 
cost of elastic deformation in the deeper layers induced 
by the buckling. The ground state geometry is thus de- 
termined by an interplay of both elastic and electronic 
effects and could be quite sensitive to different surface 
reconstructions and computational methods. 

Before we describe our own results for the ground 
state of clusters, we briefly discuss the recent literature. 
Yang et alci do not find buckling for the S19H12 clus- 
ter, either within Hartree-Fock or in a DFT calculation 
employing the B3LYP functional. Increasing the cluster 
size to Sii5Hi6, however, unambiguously shows an en- 
ergy favor for the buckled p(2 x 2) reconstruction, ~ 3-5 
kcal/mol (0.15-0.23 eV) per dimec-depending on the ba- 
sis set used. Konecny and Dorenca have used the same 
cluster and the BLYP functional to study H 2 adsorp- 
tion on Si(001)-2xl and found a buckled dimer structure 
for Si 9 Hi 2 with a = 9.6°, AE = 0.05 kcal/mol (0.002 eV) 



and d = 2.27 A while for the two-dimer cluster SiisHi6 
they obtained a = 15°, AE = 1.5 kcal/mol (0.07 eV) per 
dimer and d = 2.33 A. The two research groups use dif- 
ferent relaxation constraints and basis sets. This could 
explain the differences, given that a delicate balance of 
several effects is responsible for the ground state config- 
uration of the SigHi2 cluster. 




buckling a, (degrees) 

FIG. 2. Total energy preference AE = E(a)-E (left col- 
umn) and the gap E g {a) between HOMO and LUMO (right 
column) for the clusters as a function of dimer buckling an- 
gle a. The bottom left panel shows AE(a) for the five-layer 
p(2 x 2) slab. Eo is taken to be the energy of the corre- 
sponding unrelaxed symmetric (a — 0) cluster/slab at given 
exchange-correlation functional. Buckling is achieved by glid- 
ing the two dimer Si atoms along arcs (dashed lines) as shown 
in the inset (bottom-right panel) and keeping all other atoms 
fixed (shaded circles). 



In the present study, the ground state of the clusters 
is determined in a two-step procedure. In the first step, 
the Si dimers are tilted as a whole, while keeping the 
other cluster atoms fixed and preserving the length of 
the dimer backbonds. We sample the total energy as a 
function of the dimer buckling angle a. The results for 
LDA and the PW91 functional are shown in Fig. ||. They 
are summarized in the 'pre-relaxation' column of Table || 
for all functionals used in the calculations. 

The AE vs. a curves are mainly affected by the clus- 
ter size, rather than by the approximation used for ex- 
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TABLE I. Parameters of the clusters and the slab for dif- 
ferent exchange-correlation functionals; the energy favor for 
buckling AE is given in kcal/mol per dimer and the gap E g 
between HOMO and LUMO in eV. d refers to the Si-Si dimer 
bond length in A and a — to its buckling angle in degrees. 
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0.86 
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2.47 
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5.1 
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Sig 
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0.0 


0.96 
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0.1 


1.10 


2.24 


Sirs 
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0.5 


0.43 


15.7 


3.1 


0.88 


2.32 
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1.7 


0.33 


18.0 
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0.74 


2.40 
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0.1 


0.51 


12.2 
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0.72 


2.33 


Si 2 i 


6 


0.8 


0.27 


16.6 


3.5 


0.72 


2.44 


slab 


10 


0.8 




16.6 


4.7 




2.43 



a a and d for the Si2iH2o cluster correspond respectively to 
the buckling angle and bond length of the middle dimer. 

change and correlation. Upon increasing the cluster size, 
the minima become well pronounced and are shifted to 
larger a values, the upper bound being set by a* .In 
the case of SigHi2, a = is the only minimum when us- 
ing the BP and BLYP functionals. The values obtained 
within LDA and PW91 for this cluster are small. Thus 
the SigHi2 cluster gives no conclusive answer to the ques- 
tion whether the Si surface dimers are buckled or not. 
The larger clusters, however, clearly show a preference 
for buckling. Since we have frozen the elastic degrees of 
freedom, the only driving force that can lead to dimer 
buckling is rehybridization and charge transfer. The pre- 
relaxation study supports the view that surface dimer 
buckling is mainly driven by electronic effects, while elas- 
tic effects are responsible for the alternation of the buck- 
ling angle. 

The electronic structure of the cluster is influenced 
substantially by the dimer buckling. As a measure of 
the differences, we report the splitting between the high- 
est occupied and lowest unoccupied molecular orbitals 
(HOMO and LUMO) of the cluster, E g , in the right col- 
umn in Fig. ^J. For SigHi2 E g is a monotonically increas- 
ing function of a. It reaches its maximum for the largest 
used value of a — 20°, where the two dimer silicons and 
the two second layer neighbors of the buckled-down Si 
atom are almost coplanar. We attribute the opening of 
the HOMO-LUMO gap to changes in the symmetry char- 



acter of the orbitals. While at a = the HOMO and 
LUMO orbital correspond to the ir and tt* orbital of the 
Si dimer, they gradually develop into orbitals localized at 
one side of the Si dimer for increasing a. The orbital at 
the lower Si atom acquires more p-character and is lifted 
in energy, while the orbital at the higher Si,atom gets 
more s-character and is energetically loweredEJ. With in- 
creasing cluster size E g is strongly reduced, while the de- 
pendence on a is almost unchanged. For the larger clus- 
ters we observe that the HOMO (LUMO) wave functions 
are no longer localized at a single Si atom, but are linear 
combinations of the dangling orbitals of all the buckled- 
up (buckled-down) Si atoms of the cluster. In the Si2iH2o 
cluster, this leads to a splitting of the HOMO and LUMO 
levels into symmetric and antisymmetric states with re- 
spect to the mirror plane in the cluster Fig. |^. Thus this 
cluster reflects already to some extent the dispersion of 
surface bands observed in the slab. We expect that these 
changes in the electronic structure will also be reflected 
in the chemical reactivity, i.e. in the adsorption barrier 
£;ads f or jj 2 mo lecules. It appears that a representation of 
the surface band structure on the cluster level is required 
for a correct description of the reaction energetics. 



HOMO LUMO+ 1 




HOMO-1 _ LUMO 




FIG. 3. Top view of the SiaiHao cluster frontier orbitals. 
The contour plots are taken in the (001) plane 0.9 A above the 
buckled-up Si atom of the middle dimer (dimers are denoted 
by shaded circles). The wave functions are real- valued. Full 
(dashed) contour lines indicate positive (negative) sign. 



Relaxing the first and second silicon layers does not 
qualitatively change the results from the pre-relaxation 
study. However, the ground state energy and geome- 
try differ substantially. The buckled dimer configura- 
tions of the Sii5Hi6 and the Si2iH2o clusters are now 
more stable by 1.5-3 kcal/mol (0.07-0.13 eV) and 3.5- 
4.5 kcal/mol (0.15-0.20 eV) per dimer, respectively. The 
HOMO-LUMO gap is found to decrease with increasing 
cluster size. The reduction can be partly attributed to 
a splitting of both the HOMO and LUMO state due to 
linear combinations of dangling bonds at neighboring Si 
dimers in the larger clusters (see Fig. 3). However, a 
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decrease is also observed, in particular at a = 0, for a 
suitably averaged gap (E g ), defined as 

1 N ^ 

( E g( a )) N = JjY1 I £ lumo+» " E Houo-A a )l 
3=0 

where N is the number of dimers in the cluster. The 
reduction of the gap must be attributed to a weakening 
of the 7r-bonding with increasing cluster size. As could 
be anticipated, Si2iH2o gives closest agreement with the 
five-layer p(2 x 2) slab. The SigHi2 cluster, on the con- 
trary, is not large enough to model the properties of 
the clean Si(001)-2xl surface. Buckling in this case is 
strongly influenced by relaxation constraints and the ap- 
proximation to exchange and correlation. While for the 
two larger clusters all functionals used here give simi- 
lar results, their predictions for the SigHi2 ground state 
symmetry are qualitatively different. Though LDA and 
PW91 favor buckling with a w 7° and 5°, respectively, 
\AE\ is too small (~ 0.1 kcal/mol) to allow us to conclude 
unequivocally about the SigHi2 cluster symmetry. 

For the two-dimer Sii5Hi6 cluster using the BLYP 
functional we get d and AE values identical to those of 
Konecny and Doren, but their predicted SigHi2 geometry 
is at variance with ours. Possibly the overestimation of 
the surface relaxation due to the absence of any geomet- 
ric constraints in Ref. ^ is more crucial for the SigH^ 
cluster than for the larger ones. 

At our imposed relaxation constraints the BP and 
BLYP functionals are found to gpifi_results for SigHi 2 
in agreement with the CI methodsE3o. For the Sii5Hi 6 
cluster, however, buckling is always energetically favor- 
able within DFT, whereas recent multi-reference CI cal- 
culationsHil find the symmetric cluster to be lowest in 
energy. At present it is not clear if this is due to a lack 
of the CI calculations to recover the full correlation en- 
ergy, or due to an inadequacy of the exchange-correlation 
functionals we are using. For the largest Si2iH 2 o cluster, 
we are not aware of CI calculations addressing the dimer 
buckling. 

IV. REACTION ENERGETICS 

The energetics of dissociative adsorption and associa- 
tive desorption of H2 is characterized by three points 
along the reaction pathway, the energies of the struc- 
tures corresponding to the monohydride phase En, the 
transition state energy _&rs an d the sum of bare cluster 
energy Eoq and that of the free H2 molecule. Hence, the 
quantities of interest are defined by the differences 

gads _ _ _j_ 22(H 2 )), 

tj a — ^TS — 

-Erxri = -EoO + -E , (H2) — E\\. 




0.8 1 1.2 1.4 

H1-H2 distance (A) 



FIG. 4. Sketch of the ridge method search for TS and the 
adiabatic PES for the Sig cluster within LDA. Crossed circles 
(©) denote the successive approximations to the TS. Contour 
spacing is 0.05 eV and energy is measured with respect to 
Soo + -E(H2) (solid contour line). Transition state geometry is 
shown in the inset. PES is plotted for H2 molecule impinging 
perpendicularly to the surface with its center of mass right 
above the Sil site, Z being the distance between them. 

The monohydride geometry of each cluster is obtained 
by saturating the dangling bonds of one Si dimer with 
H atoms. The equilibrium Si-H bond is 1.52 A within 
LDA, PW91 and BLYP. When using the BP functional 
a bond length larger by 0.03 A is obtained. The changes 
in d of the hydrogenated dimer are mainly governed by 
the approximation to the exchange and correlation, while 
the cluster size accounts only for differences of - 0.01 A 
for a given functional. We note that our TS and mono- 
hydride structures imply reaction of a single H 2 molecule 
with the Si(001) surface. For comparison to experimen- 
tal data taken at finite hydrogen coverage, information 
about the coverage dependence of the energetics is also 
required. A single H2 molecule per p(2 x 2) unit cell in the 
extended slab case corresponds to a coverage of 8 = 0.5 
monolayers (ML). To assess the coverage dependence of 
the reaction energy E rxn , test calculations were carried 
out for a completely covered slab, 8 = 1 ML. To obtain 
a better understanding of the role of buckling for the 
reaction energetics, we also performed test calculations 
for the symmetric p(2 x 1) reconstructed Si(001) surface 
as possible reference of the clean surface. We note that 
the occupation of adjacent Si dimers by hydrogens could 
have a (probably small) influence on the transition state 
geometries and energies via interaction with the neigh- 
boring monohydride. This effect can be studied to some 
extent with the help of the two- and three-dimer clus- 
ters, which allow for adjacent doubly occupied Si dimers. 



G 



TABLE II. Transition state geometry parameters in A (see 
the inset in Fig. ^). «' denotes the buckling angle in degrees 
of the unoccupied dimer(s) in the case of slab and Siis, S121 
clusters. 



p 

A H1-H2 


^Sil-Hl 




- rL Sil-H2 


D 

XX Si2-H2 


d 


c* TS 


a' 


LDA 
















Sig 


0.88 


2.12 


2.10 


2.38 


2.30 


13.7 


— 


Siis 


0.97 


1.93 


1.98 


2.24 


2.38 


12.1 


15.6 


Sin. 


1.06 


1.95 


2.14 


2.09 


2.38 


13.0 


15.5 


slab 


0.93 


2.00 


2.11 


2.24 


2.39 


13.5 


19.0 


BP 
















Si 9 


0.94 


1.96 


1.95 


2.34 


2.42 


12.2 


— 


Siis 


1.02 


1.84 


1.95 


2.19 


2.46 


11.5 


15.5 


Si 2 i 


1.02 


1.82 


2.01 


2.12 


2.47 


11.6 


14.8 


slab 


1.03 


1.80 


1.94 


2.11 


2.50 


13.1 


18.7 


PW91 
















Sig 


0.89 


1.98 


1.96 


2.32 


2.37 


12.2 




Siis 


0.99 


1.84 


1.94 


2.19 


2.40 


10.7 


14.6 


Siai 


1.08 


1.89 


2.11 


2.04 


2.41 


13.6 


14.9 


slab 


0.99 


1.83 


1.99 


2.12 


2.43 


12.7 


18.4 


BLYP 
















Sig 


0.91 


1.93 


1.95 


2.31 


2.40 


11.9 




Siis 


0.96 


1.86 


1.98 


2.23 


2.44 


10.2 


12.6 


Siai 


0.96 


1.83 


2.02 


2.11 


2.48 


12.1 


13.8 


slab 


1.02 


1.74 


1.89 


2.12 


2.48 


11.1 


18.4 



To get some insight into the influence of neighboring 
monohydrides on the asymmetric TS, we have performed 
calculations with SiisHx6+ x and Si2iH 2 o+a; clusters with 
x = 4 and 6, respectively. 

In this study, we concentrate on the pre-pairing sce- 
nario for the H 2 reaction with the Si(001) surface. Con- 
sequently, we locate the asymmetric TS of H2 desorption 
from a single Si dimer for all clusters and functionals 
used in this study. In principle this can be achieved by 
mapping out the related potential-energy surface (PES), 
like e.g. in Ref. [l5|. As an example, the potential energy 
as a function of the distance between the two H atoms 
and the HVcluster distance Z is shown in Fig. ^. For 
each configuration of the two H atoms, the Si atoms in 
the two topmost layers have been relaxed. They follow 
'adiabatically' the motion of the H atoms. Thus we make 
sure that the lowest possible TS in the multidimensional 
space of all mobile atomic coordinates is found. An alter- 
native and generally faster approach uses a search algo- 
rithm, thus avoiding the need to map out all the points 
in the PES. All degrees of freedom of the H2 molecule 
and the topmost two cluster/slab Si. layers are included 
in the search. The ridge mcthodEj implemented here 
starts from a pair of coordinates xq and x\ in the multi- 
dimensional configuration space of the system, which de- 
note the reactants (H2 above the bare surface) and the 
products (the monohydride), respectively. An iterative 
search is then performed by halving the interval [xq,x5J 
in each step. To reduce the number of required steps, we 
have shifted the input coordinates towards some initial 
guess (xq, x\) closer to the TS. A projection of the search 



path onto a two-dimensional slice of the coordinate space 
is shown in Fig. |i| for the S19H42+2 cluster within LDA. 
As seen, there is excellent agreement between the two 
schemes. 

The geometries of the calculated TS are collected 
in Table |l[ As expected, the H-H bond is stretched in 
the asymmetric TS configuration of all Si x Hj, + 2 clusters, 
with -R H1 _ H2 being largest for dissociation on the Si 2 i 
cluster. As a consequence, the atom H2 at the transition 
state is by 0.2-0.3 A closer to the buckled-up Si2 atom 
as compared with the Sig cluster. The presence of two 
H atoms over the Sil site partially blocks the mechanism 
that leads to anticorrelated dimer buckling and therefore 
ev TS < a. For the Si2i cluster the unoccupied dimers 
are somewhat affected by the adsorption event, but their 
buckling angle a! is only a few degrees smaller than that 
of the clean surface. This small change is due to the 
dimer row termination by hydrogens in the clusters, and 
is absent in the slab geometries, where a' « a. 




Si 15 Si 21 
size of the system 



slab 



FIG. 5. H2/Si(001) reaction energetics as a function of the 
size of the system used to model the Si(001)-2x 1 surface (see 
also Table m). 



Most noticeable in the results shown in Fig. || (nu- 
merical values are compiled in Table [II) is a clear de- 
pendence of E^ ds , Eq QS and E rxn on cluster size. The 
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TABLE III. Activation and reaction energies in kcal/mol 
(eV) per molecule for H2 adsorption/desorption on 
Si(001)-2xl surface via pre-pairing mechanism (ZPE not in- 
cluded) . 









^rxn 


LDA 










Sig 


10.2 (0.44) 


64.9 


(2.81) 


54.7 (2.37) 


Siis 


5.7 (0.25) 


58.4 


(2.53) 


52.7 (2.29) 


Siai 


3.5 (0.15) 


53.1 


(2.30) 


49.6 (2.15) 


slab 


0.9 (0.04) 


49.0 


(2.12) 


48.1 (2.09) 


BP 










Si 9 


18.1 (0.78) 


67.3 


(2.92) 


49.3 (2.14) 


Siis 


13.3 (0.56) 


61.4 


(2.66) 


48.1 (2.09) 


Si2i 


10.2 (0.44) 


53.9 


(2.34) 


43.8 (1.90) 


slab 


11.6 (0.50) 


53.4 


(2.32) 


41.8 (1.81) 


PW91 












1 k n ( r\ an\ 
lo.y (U.DyJ 


66.0 


(2.86) 


OU.l (A. LI) 


Siis 


14.8 (0.64) 


63.0 


(2.73) 


48.2 (2.09) 


Siai 


12.0 (0.52) 


56.7 


(2.46) 


44.7 (1.94) 


slab 


9.6 (0.42) 


52.5 


(2.28) 


42.9 (1.86) 


BLYP 










Sig 


20.2 (0.88) 


72.8 


(3.16) 


49.7 (2.16) 


Siis 


15.8 (0.69) 


64.6 


(2.80) 


48.8 (2.12) 


Si2i 


14.5 (0.63) 


58.8 


(2.55) 


44.3 (1.92) 


slab 


15.5 (0.67) 


58.4 


(2.53) 


42.9 (1.86) 


common 


trend for all functionals 


is an 


effective flatten- 



ing of PES along the reaction path with increasing cluster 
size. Both the reaction energy and the adsorption barrier 
are reduced. In general, for all clusters considered here, 
LDA gives a lower bound for the activation barriers and 
an upper one for i? rxn . For a given size of the cluster, 
we can compare the effect of the exchange-correlation 
functional employed in the calculation. As far as E mi is 
concerned, all gradient-corrected functionals behave in a 
very similar way (see bottom panel of Fig. ^|). This gives 
credibility to the statement that gradient-corrected DFT 
yields an accurate description of reaction energies. The 
quantities E^ ds , E dcs , involving transition state energies, 
show a stronger variation. The differences between BP 
and PW91 are only ~ 2 kcal/mol, with the sign depend- 
ing on the size of the system. The BLYP functional gives 
the highest value for E^ ds ' dca of all tested functionals. 
Since similar performance has been already established 
for small systemaid and the Sig clustered, our calculations 
confirm these results for extended clusters and slabs. 

The Si2i cluster displays both the lowest adsorption 
and desorption barriers. Comparison between this clus- 
ter and the slab shows that their predictions agree to 
within ~ 3 kcal/mol for all quantities. We conclude that 
the S121 cluster gives a fair description of the Si(001) sur- 
face, while the others are inadequate approximations for 
the surface. While the desorption barriers derived from 
the S121 cluster are in the range of 56 ± 3 kcal/mol for 
all functionals, E^ ds is more sensitive to the functional 
used, with values covering a range of 11 kcal/mol. The 
adsorption barriers derived from the gradient-corrected 



functionals using this cluster are 7-11 kcal/mol higher 
than the LDA barrier. This is in accord with the estab- 
lished picture that LDA— tends to underestimate adsorp- 
tion barriers at surfaces^. We note that the S121 cluster- 
yields higher barriers for the PW91 functional than for 
BP, while the slab calculations give the opposite result. 




(a) (b) 

FIG. 6. Total valence electron density n(r) (grey shading) 
and the density difference An(r) (contour plot) in the plane 
containing H2 molecule and the Si-Si surface dimer at the 
TS geometries of the Si 9 Hi2+2 (a) and the Si2iH2o+2 cluster 
(b) calculated with the PW91 functional. Density difference is 
defined by An(r) = n(r) -n^ st (r) -n^|(r). The full contour 
lines correspond to An > and dashed lines to An < 0. The 
plot levels are the same for both clusters. 



It is instructive to analyze why the smaller Sig and 
Siis clusters give a poor description of the physics at 
the Si(001) surface, despite the fact that the Si-H bound 
is localized and should therefore be well represented al- 
ready in the smallest cluster used. For E rxn , the variation 
with cluster size is mainly due to differences in the bare 
clusters which are used as reference states. The buck- 
ling of the Si dimers characteristic for the Si(001) sur- 
face only fully develops in the largest clusters, because 
it requires a correct description of the electronic surface 
states. Thus non-local electronic effects enter the calcu- 
lation of the reaction energy. They give rise to a reduc- 
tion of -Erxn by about 5 kcal/mol compared to the value 
obtained with the Sig cluster. For E^ da , the differences 
between small clusters and the slab calculation are even 
larger. This is due to the fact that the dangling bonds 
of the Si dimer act as frontier orbitals in the reaction 
with H2. The importance of the change in the HOMO 
and LUMO position that accompany dimer buckling has 
also ibf en emphasized-in recent first principles-studies of 
H 2 OEl, C2H2N, BH3EI and 1,3-cyclohexadienEl reactions 
with the Si(001) surface. As our calculations show, the 
HOMO-LUMO gap E g is substantially reduced when go- 
ing from the (mainly 7r-bonded) Si dimer in Sig to the 
surface band states in the slab. At the transition state 
for adsorption, the H2 molecular orbitals, in particular 
the antibonding 3 E U orbital, mix with both the HOMO 
and LUMO orbital, which is facilitated by a small E g . At 
the TS of the Sig cluster, the H 2 mainly interacts with 
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TABLE IV. Coverage dependence of the reaction energy 
E (9) for the slab. 







(9), kcal/mol (eV) per H 2 


0.5 ML 


1 ML 


1 -> 0.5 ML 


LDA 


48.1 (2.09) 


49.7 (2.16) 


51.3 (2.22) 


BP 


41.8 (1.81) 


43.8 (1.90) 


45.7 (1.98) 


PW91 


42.9 (1.86) 


44.7 (1.94) 


46.5 (2.02) 


BLYP 


42.9 (1.86) 


44.3 (1.92) 


45.7 (1.98) 



one dangling bond of the lower Si atom (see Fig. In 
the more extended systems with smaller E g , there is also 
a direct interaction between the H2 molecule and the up- 
per Si atom. This shows up in the more localized induced 
charge distribution between H2 and Si2 (right panel of 
Fig. [j|), and in the different geometrical structure of the 
transition state. 

The trends in the reaction energetics outlined above 
do not change substantially if TS and monohydride con- 
figurations of the Sii5, S121 clusters with more than 
one monohydride are employed. Test calculations using 
Sii5Hi6+4 and Si2iH2o+6 clusters show that finite cover- 
age effects at the TS introduce a variation of less than 
3 kcal/mol in the calculated barriers. The-coverage de- 
pendence of .Erxn was studied for the slabH with initial 
coverage = 1 ML, Table IV. It is energetically more 
expensive to desorb a H2 from one of the dimers if the 
other stays monohydride. Thus, E rxn (l — * 0.5 ML) is 
about 3-4 kcal/mol larger than E rKn calculated for ini- 
tial coverage 6 = 0.5 ML. The reaction endothcrmic- 



ity for the removal of a whole monolayer comes out as 
the average of the endothermicities associated with the 
removal of each of the two monohydrides in the unit 
cell. It is interesting to note that E lxn for the Si2i clus- 
ter with one monohydride (see Table III) coincides with 
E rxn (l ML), rather than with the low-coverage limit of 
the slab, l? rxn (0.5 ML). The coverage dependence cal- 
culated for the slab shows the same trend as the ex- 
perimental data by Flowers et a/O, where a slight in- 
crease of the desorption energy with coverage is observed, 
E^ CS (Q) = (55.8 + 1.1 x 9) kcal/mol for initial coverages 
in the range 0.01-1 ML. 

The reaction energy E rxn calculated with respect to 
the symmetric 2x1 reconstructed Si(001) surface shows 
very similar behavior. The respective E rxn (Q) values are 
increased by a few kcal/mol. In this case, when a H2 
molecule desorbs, it leaves behind a symmetric unoccu- 
pied dimer and therefore the final state is ~ AE higher 
in energy than for the buckled surface reconstruction. 



the cluster predictions. All calculations were performed 
within LDA and with the non-local BP, PW91 and BLYP 
exchange-correlation functionals. As our results show, a 
conservative conclusion can be drawn that the most fre- 
quently used SigHi2 cluster is not large enough cither 
to model the properties of the bare Si(001) surface or 
the molecular H2 dissociation on it. The latter stems 
from the fact that this cluster is not capable of recov- 
ering the surface electronic structure. Though the H2 
reaction with the Si(001) surface is considered to be a 
highly localized event, non-local effects enter the reac- 
tion energetics via their influence on the surface bands. 
Hence, one could also expect different performance for 
the various functionals. 

Our analysis shows that the quality of a given 
exchange-correlation functional should not be assessed 
without referring to the particular size of the cluster em- 
ployed to study the H2/Si(001) adsorption/desorption 
process. Indeed, it is evident, by inspecting the BLYP 
section of Table III for example, that one could infer for 



the Sig cluster an activation barrier to desorption much 
higher than the experimental values. Hence, as usually 
proceeded, the pre-pairing mechanism could be ruled out 
on energetic grounds. Such a conclusion, however, seems 
to be premature if one refers to the slab or, eventually, 
the Siai cluster prediction within the same functional. 

In contrast to the one-dimer cluster approximation, the 
Si2iH2o cluster was found to be close in its predictions to 
the slab calculations in all aspects of the reaction consid- 
ered. With E^ cs — 56 ± 3 kcal/mol it is also well within 
the range of the experimentally determined desorption 
barriers. The relatively large spread in the calculated en- 
ergies comes from the functionals used, with the BLYP 
functional giving the largest value for £:ads,dcs^ ZPE cor- 
rections were not considered, but as they amount to a 
few kcal/mol at most and are essentially the same for 
different functionals, no qualitative change of our results 
is expected upon their inclusion. Concluding, our find- 
ings suggest that the effect of the cluster size in modeling 
the H2 reaction with the Si(001) surface is significant and 
some of the previous works may well need a revision. 
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V. CONCLUSIONS 



We have presented a systematic ab initio study of the 
H2/Si(001) reaction energetics employing three clusters 
in plane-wave DFT calculations. A five-layer p(2 x 2) 
slab was used as a reference to analyze the convergence of 
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